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ABSTRACT 

We study the conditions for disk galaxies to produce superbubbles that can break out 
of the disk and produce a galactic wind. We argue that the threshold surface density of 
supernovae rate for seeding a wind depends on the ability of superbubble energetics to 
compensate for radiative cooling. We first adapt Kompancets formalism for expanding 
bubbles in a stratified medium to the case of continuous energy injection and include the 
effects of radiative cooling in the shell. With the help of hydrodynamic simulations, we 
then study the evolution of superbubbles evolving in stratified disks with typical disk 
parameters. We identify two crucial energy injection rates that differ in their effects, 
the corresponding breakout ranging from being gentle to a vigorous one. (a) Superbub- 
bles that break out of the disk with a Mach number of order 2-3 correspond to an 
energy injection rate of order 10~* erg cm~^ s~^, which is relevant for disk galaxies 
with synchrotron emitting gas in the extra-planar regions, (b) A larger energy injection 
threshold, of order 10"'^ erg cm^^ s^^, or equivalently, a star formation surface density 
of ^ 0.1 Mq yr~^ kpc~^, corresponds to superbubbles with a Mach number ~ 5-10. 
While the milder superbubbles can be produced by large OB associations, the latter kind 
requires super-starclusters. These derived conditions compare well with observations of 
disk galaxies with winds and the existence of multiphase halo gas. Furthermore, we find 
that contrary to the general belief that superbubbles fragment through Rayleigh- Taylor 
(RT) instability when they reach a vertical height of order the scale height, the super- 
bubbles arc first affected by thermal instability for typical disk parameters and that the 
RT instability takes over when the shells reach a distance of approximately twice the 
scale height. 
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1 INTRODUCTION 

Observations of nearby and high-redshift galaxies have shown 
that star formation in them often leads to galactic winds. 
Starburst galaxies, with star formation rate (SFR) in excess 
of a few tens of Mq yr~^ are known to excite such outflows. 
However, Heckman (2002) pointed out that it is not the av- 
erage SFR, but the SFR surface density which is a deciding 
factor for the existence of outflows. He found a threshold SFR 
surface density of ~ 0.1 M0 kpc~^ yr~^ as a prerequisite for 
starbursts to be able to produce galactic winds. 
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The standard scenario of star formation leading to the 
wind phenomena posits that super-starclusters give rise to a 
large number of supenovae (SN) in a relative small region, 
which can produce a superbubble in the disk and can break 
out of the disk with enough momentum to produce a wind. 
Such super-star clusters, or young globular clusters, have been 
observed to have masses in the range of few X 10^-6 X 10^ 
M0 within a typical radius of ~ 3-10 pc (Ho 1997; Marti'n- 
Hernandez et al. 2005, Walcher et al. 2005). The large amount 
of energy deposited into the interstellar medium (ISM) by 
these objects in the form of UV radiation and mechanical en- 
ergy is believed to be an important feedback process. The me- 
chanical energy from these super-starclusters has been shown 
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to be important for the superbubble produced by the com- 
bined SNe to break out of the disk and produce a large scale 
wind (e.g., Tenorio-Tagle, Silich, Mufioz-Tunon 2003). 

There have been a number of calculations, both ana- 
lytical and numerical, dealing with the breakout of super- 
bubbles from disk galaxies. The conditions for breakout de- 
pend strongly on the assumption of the stratification of gas 
in the disk. Consider an exponentially stratified disk with 
mid-plane ambient gas pressure Pq, gas density po, scale 
height zo, and a bubble being blown by mechanical lumi- 
nosity C Mac Low & McCray (1988) defined a dimensionless 
parameter D = Cp^^^ / {Pq^^ Zq) , and noticed in their nu- 
merical simulations that breakout of bubbles occurred when 
D ^ 100. The importance of this parameter can be under- 
stood by considering the self-similar evolution of a super- 
bubble driven by an energy injection rate of £, given by 
r ~ (£^Vpo)'''^ and r ~ (3/5)(£/po)'''^r This im- 
plies a speed of ~ {3/5){£/ pqZq)^^^ oc D^^^ when the su- 
perbubble reaches a distance of the scale height, for an am- 
bient gas at a given temperature. According to this criterion, 
for a scale height zo = 200 pc, and mid-plane gas density 
p.Hno ~ 2.3 X 10"^" g cm"^, Po/h ~ nolO* K cm"^, a bub- 
ble with total mechanical luminosity oi C ^ 3.8 x 10"^^ erg s^^ 
will be able to breakout of the ISM. 

Basu et al. (1999) defined a dimensionless parameter 
b = {27/154Tvy/'^/:^^^pl^*Po'^^^Zg^ which is a ratio of the 
radius where the Mach number of the superbubble becomes 
unity, to the scale height. This is motivated by the self-similar 
solution of a stellar wind, r ~ {125/15471^^"^ p^^^^t^^^ . They 
showed that this parameter is related to the above mentioned 
D parameter as D = 17. 9&^. In other words, a superbubble 
with & < 1 is likely to be confined where as blowout will occur 
for 6^1. 

Koo & McKee (1992) analytically determined a con- 
dition for the breakout. Since the bubbles accelerate after 
reaching a distance of order the scale height, owing to the 
rapidly decreasing density, it becomes liable to fragment due 
to Rayleigh- Taylor instability. If the Mach number of the bub- 
ble at scale height is > 3, then they argued that the bubble 
would be able to breakout. They used radiative bubble model 
of Weaver et al (1977) for a uniform density atmosphere in 
order to derive a critical mechanical luminosity for which the 
Mach number is unity, Ccr ~ 17.9poZqcI, where Cs is the 
isothermal sound speed of the ambient gas. The Mac Low 
& McCray condition of D ^ 100 translates to £/Ccr ^ 5. 
As we will find in our simulations, the Mach number of a 
bubble after breakout is of order (l/5cs){£/poZo)^^^ . There- 
fore, the Mac Low-McCray condition of _D ^ 100 translates 
to the condition that the Mach number at breakout is of or- 
der unity. We also note that they considered superbubbles 
that originated at a height from the mid-plane, which made 
it easier for bubbles to break out. Our simulations show that 
the critical luminosity for Mach number at a distance of the 
scale height to be unity is Ccr ^ 125poZqC'1, larger than the 
estimate of Koo & McKee (1992). 

Koo & McKee (1992) then considered an additional 
strata of HII gas with a scale height of 1 kpc and mid-plane 



number density 0.025 cm ^, and found the breakout condition 
to be of order Nqb ^ 800, or equivalently, £ ^ 4.14 x 10^* 
erg s~^. Silich & Tenorio-Tagie (2001) considered the effect of 
halo gas pressure and determined a minimum energy for the 
superbubble to blow out of the galaxies (with both disk and 
spherical ISM distribution) with ISM gas mass in the range of 
lO'^-lO^ Mq. For a disk galaxy with Mism ~ 10^ Mq, they 
found a minimum energy of ~ lO'** erg s~^, corresponding to 
iVoB ~ 1000. 

As Heckman (2002) has emphasized, it is the surface den- 
sity of SFR that determines the condition for the existence 
of galactic winds, and not the total luminosity. To translate 
the above energy conditions into a surface density, we need 
to estimate the surface area of such bubbles at the breakout 
epoch. In this paper, we re-visit this issue in order to under- 
stand the empirical threshold SFR surface density for galactic 
winds. Murray, Menard, Thompson (2011) have recently ar- 
gued that radiation pressure from UV radiation from a disk 
with a SFR surface density larger than 0.1 M© kpc~^ yr~^ 
can produce a large scale wind. This estimate however cru- 
cially depends on the assumption of the grain opacity, and as 
Sharma & Nath (2012) have shown the relevant opacity at 
UV may fall short of the requirements. 

There have also been studies on the existence of multi- 
phase gas in the halos of spiral galaxies, and their connec- 
tion to the star formation properties in the disk. Dahlem, 
Lisenfeld, GoUa (1995) considered nine edge-on galaxies with 
extended synchrotron emitting halo gas, and derived a min- 
imum value of surface density of energy injection for super- 
bubble breakout, as ~ 10^** erg s^^ cm^^. Tiillmann et al. 
(2006) further considered X-ray, radio and far-infrared (FIR) 
emission from the extended halo gas in a sample of 23 edge-on 
spiral galaxies, and found that the halo contained gas at low 
and high temperatures (multiphase) if the surface density of 
energy injection in the disk exceeds ~ 10""^ erg s~^ cm~^. 
If the existence of multiphase halo gas depends on the pro- 
cess of superbubbles breaking out of the disk and depositing 
hot interior gas ( as suggested by Tomisaka & Ikeuchi 1986; 
Tenorio-Tagle, Rozyczka, Bodenheimer 1990), as well as cold 
gas in the shell, then it would be interesting to compare the 
energetics of such superbubbles and the observed threshold 
energy injection rate. 

In this paper, we study the standard scenario of thermal 
pressure of the gas interior to superbubbles being the driving 
mechanism for the wind, and derive a threshold condition 
for the superwind. We find that radiative loss of energy is 
important for the dynamics of shocks, and the inclusion of 
radiation loss increases the energy budget for the bubbles to 
breakout of the disk and produce a wind. We also find that our 
estimate of the threshold energy requirement can explain the 
observed threshold SFR surface density for galactic outfiows. 

The paper is organized as follows. In §2 we derive an 
order-of-magnitude estimate of the threshold based on the 
key idea that the superbubble energetics needs to balance 
radiative cooling. Then we present the analytical formalism 
in §3 and discuss the results in §4. We then present the results 
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from numerical simulations in §5, and discuss the effect of 
thermal and RT insability in §6. 



2 ANALYTIC ESTIMATES 

To begin with, we derive a threshold rate of SNe for a super- 
bubble to continue to grow and ultimately breakout of the 
disk from simple arguments. We can first consider the condi- 
tion that the superbubble is able to drive a strong shock in 
the disk. This requires the volume energy injection time scale 
to be shorter than the sound crossing time. In other words, 
if we consider a region of radius R in the disk and an energy 
injection rate of £, then one needs 



< R/c, , 



(1) 



l.5nkT 
C/iA-RR^S) 

where Cs is the sound speed. This gives a lower limit of 
C/{itR^) > 3 X IQ-'^nTl''^ erg s"^ cm"^, where n is the 
ambient gas particle density in cm" "^ and T = r4 lO" K. 

A second, and more stringent, constraint on SNe luminos- 
ity comes from accounting for radiative losses. Let us assume 
that when a SN remnant enters the radiative stage it quickly 
looses its energy and does not contribute to the energy in- 
put of the superbubble. Assume then that the radiative stage 
begins when the post-shock temperature is Ts — 10^ K such 
that the radiation loss function is maximum and much larger 
than the minimum at ~ 10^ K. We therefore define the time 
when a SN remnant losses its energy as time when the shock 
velocity is Vg = 70 km s~^ (corresponding to the post-shock 
temperature of 10^ K). It determines the corresponding time 
and radius as 

p,l/3 p,l/3 

ta = 3xl0"^yr, i?a = 50^pc. (2) 
One can therefore define the coherency condition as, 

— RitaVs^ > 1 , (3) 

which means that before a SN remnant stalls because of cool- 
ing losses, another SN explosion injects energy into the rem- 
nant and forms a single bubble. This condition determines 
the required SN rate 



> 6 X 10 



4/3 

-12 / " \ o^T -1 -3 

— SNe yr pc 

-C/51 . 



(4) 



We can estimate the surface density of SNe, by mul- 
tiplying this rate density by the scale height, which is 
the height of a bubble at the epoch of breakout. For a 
scale height of 500 zq = 2:0,0.5 pc, this corresponds to 
3 X 10-^(n/£;5i)*''^ 20,0.5 SNe yr"^ kpc'^ (The scale height 
is relevant here because, as we will see later, the maxi- 
mum radius of bubbles in the plane parallel to the disk 
is of order nzo.) Finally, we recall that for a Salpeter 
IMF, one SN corresponds to 150 M0 of stellar mass, con- 
sidering stars in the range of 1-100 M©. Therefore the 
threshold condition for SFR surface density becomes ~ 



0.5(n/S5i)*''3 



2o,o.5M0 yr 



^ kpc ^. The corresponding sur- 



face density of energy injection is ~ 0.01 n^^^ 20,0.5 
erg s"^ cm~^. It is interesting to find that these above esti- 
mates of the threshold energy injection or SFR surface density 
are comparable to the observed threshold for the existence of 
multiphase halo gas (Tiillmann et al. 2006) and superwinds 
(Heckman 2002). 



3 KOMPANEETS APPROXIMATION 

We first discuss the expansion of blastwaves in a stratified at- 
mosphere, in the adiabatic case and then for radiative shocks. 
Kompaneets (1960) had first analytically worked out the case 
of adiabatic shocks in this case (see, e.g., Bisnovatyi-Kogan 
& Silich 1995). Consider an exponentially stratified medium 
described by p{z) = po exp ( — 2/20), where po is the midplane 
density and 20 is the scale height and Eq is the explosion en- 
ergy. It is assumed that the shock pressure is uniform, and is 
given by. 



(7 - l)\Eo 
V 



(5) 



where A ~ 1 (Kompaneets 1960) is a constant that differenti- 
ates the shock pressure from the average pressure inside the 
bubble; Bisnovatyi-Kogan & Silich (1995) evaluated A = 1.33 
which we adopt here. We define a dimensionless time-like pa- 
rameter as. 



y ■- 



(7" - l)A£:o 
2poV 



dt. 



(6) 



Where Eth is the thermal energy of the interior gas, V is the 
volume of the bubble and t is the time. The shape of the shock 
front is derived as, 



r = 220 arccos-l i exp (2/220) 



y 



-f exp (-2/20) 



(7) 



The location of the top and bottom of the bubble then follows 
by setting r — Q [ with y = y/zo), 



z±{y) = -22o.ln(lTy/2), 



(8) 



which shows that the top of the bubble reaches infinity when 
y — > 220 while t remains finite. This implies that the bub- 
ble accelerates in the 2-direction due to stratification, after 
an initial deceleration phase when the bubble is small and 
spherical, as in the usual Sedov- Taylor solution. The maxi- 
mum cylindrical radius of the bubble is also obtained from 
the above solution by putting (dr/dz) = 0, 



'"max(y) = 220 arcsin {y/2) . 



(9) 



The 2-component of the velocity of the topmost point of the 
bubble is given by. 



(1 - y/2) 



(72-1) A£;o 



PoV{t) 



(10) 
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Figure 1. The ratio of cooling time to time (tcool/*) is plotted 
against the height of adiabatic superbubble with continuous energy 
injection, for Nqb = 20, no = 0.1 cm~^, zq = 200 pc. 



3.1 Continuous energy injection 

We can extend Kompaneets approximation and radiative 
blastwave calculation to the case of continuous energy injec- 
tion. Schiano (1985) had done a similar calculation in the 
case of an active galactic nucleus. Consider an association 
with Nob stars with masses above 8 M0, which ultimately 
go supernovae. If we consider the main-sequence lifetime as 
rsAT ~ 5 X 10^ yr for these stars, then the total mechanical 
luminosity of the SN in the association can written as. 



£ = 6.3 X W^^'NoB Esi (tsn/S x loVO"' 



erg s 



(11) 



where supernova energy is 10^^ £51 erg. As McCray and 
Kafatos (1987) have argued, since the main sequence life time 
scales with the stellar mass as r oc since the initial 

mass function (IMF) is given by, dAf,/d(logM*) oc M*"^'^^, 
for a Salpeter IMF, the rate of SN will scale with time 

as oc oc ^2-35/1.6 ^-1/1.6^1 ^ ^1.35/1.6-1 ^^-^^ 

dA/* at ' 

roughly constant in time. Here we have used oc t~^^^'^~^, 
given the above mentioned dependence of stellar main se- 
quence lifetime. Therefore we can write, for the adiabatic case, 
the total energy in the superbubble as Eth ~ jCt. 

Instead of egn llOl the 2— velocity of the top of the bubble 
is then given by. 



1 (7^-1) xa 



(l-y/2)V 2 poV{t) 



(12) 



and the corresponding y parameter is also written in terms 
of t, as 



/ (7i-i)A£f' 



(13) 



These equations can determine the dynamics of the super- 
bubble in the case of continuous energy injection. 



3.2 Radiative loss with continuous injection 

Radiative losses can be important for the dynamics of both 
the blastwave and a superbubble with continuous energy 
injection. Shocks become radiative when the cooling time 
icooi ^ t. The cooling time behind the shell can be esti- 
mated as tcooi = 1.5fcT/(4nA(r)), for a strong shock with 
n — no exp(— z/zo), and the shock temperature being esti- 
mated from the shock speed (in the 2;— direction, say). We 
assume a cooling function, as given by Eqn 12 in Sharma et 
al. 2010, appropriate for gas with 10* ^ T{= lO'^Te) ^ 10^ 
K. Figure [1] shows the ratio t^ooi/t as a function of the bub- 
ble height z+ for bubbles with continuous energy injection 
(thick solid line), for Nqb = 20, no = 0.1 cm"^, zo = 200 pc. 
The curves shows that the shock enters the radiative phase 
much before reaching the scale height even for a low ambient 
density of no = 0.1 cm""^. 

Radiation loss from the shocked medium can therefore 
be important. Kovalenko & Shchekinov (1985) had calculated 
the dynamics of a blastwave with radiative loss, assuming that 
the shock kinetic energy is converted into thermal energy of 
gas in a thin shell behind it, and that radiative loss from this 
shell keeps the shock isothermal. It can then be shown that 
for a strong shock the energy lost per unit mass is ~ (1/2)Ms, 
where Ua is the shock speed. From the Hugoniot condition for 
a strong shock that. 



(7 + l)Ps _ (7' - l)ASt, 



2po 



2poV(t) 



(14) 



where Eth is the thermal energy of the shocked gas. The struc- 
ture of the shock in this case is such that the interior gas 
remains hot and adiabatic, whereas the shocked ambient gas 
that is swept into a shell loses its energy radiatively and is 
kept at a constant temperature ( at ~ 10* K). We note that 
Mac Low & McCray (1988) showed that the radiative loss 
from the interior hot gas of the bubble does not change the 
dynamics of the bubble. 

Following the calculation of Kovalenko & Shchekinov 
(1985) for a radiative blastwave, we assume that bubbles 
with continuous energy injection also forms an isothermal 
thin shell, after a certain time ti when it enters the radia- 
tive phase. For simplicity, we also assume a self-similar solu- 
tion for a spherical shock, of the type given by Weaver et al. 
(1977), Ta = AC^^'H'-^^'\ where A is a constant depending on 
the ambient density. 

Using the result derived in eqn [14] that the amount 
of energy lost per unit volume is (l/2)po^'s = il'^ ~ 
l)XEth/{'lV{t)), we can write for the evolution of thermal 
energy in this case. 



rr 

Eth{r) = /:f-7r(7^ - 1) / 

J ri 



Eth{r) 2 

r dr , 



V 



^2/3 / \5/3 / 2 , 

C ' (-^) ' -7r(7 - 1) 



Ethir) 
V 



r dr . (15) 



We can explicitly solve this equation for a spherical 
shock, and then use the results to estimate the z— velocity 
of an oval shaped bubble. For a spherical shock (with volume 
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10" 
time (Myr) 

Figure 2. The evolution of the ratio of to Cs (the sound speed 
for an ambient gas at 10* K) is plotted against time, for an adia- 
batic blastwave (thick solid line), adiabatic superbubble with con- 
tinuous energy injection (dashed) and with radiative loss (dotted 
line). 



V = ^Tvr^), the energy equation (no. llSp can be shown to 
yield a solution of the type Ethir) = br" , where 



fe[a+ -(7 



1)] 



r = — 



5 £2/^ 



.2/3 



(16) 



Comparing the powers of r from both the sides we get, a = |. 
Putting this value of a in eqn [TS] and comparing the coeffi- 
cients of time on both sides we get, 



9^5/3 



Therefore Ethir) becomes, 
5 



Eth{r) 



9 



-Ct, 



(17) 



(18) 



showing that roughly half of the total energy is radiated away 
and half is stored as the shell kinetic energy, compared to the 
adiabatic case in which we had assumed Eth = This is 
an asymptotic value of the loss in the limit r 2> ri, in the 
regime where the approximation _B oc is valid. We can 
therefore use equations [T^] and 1131 with the above value of 
Eth, and determine the dynamics of a radiative superbubble 
with continuous energy injection. 



4 ANALYTIC RESULTS 

Figure 2 shows the evolution of the Mach number for a 10* 
K gas as a function of time, for an adiabatic blastwave, 
a superbubble with continuous energy injection with and 
without radiative loss. It is convenient to define a dynami- 
cal time scale for this problem (Mac Low & McCray 1988), 



td 



(po/C) ' , which is the expected time to reach the 



scale height for a self-similar evolution of superbubbles. For 
zo = 200pc, £ ~ 1.3 X lO^'^ erg s"^ and po ~ lO"^''^ g cm"^ 
(for p. ~ 0.6) , td ~ 2.84 Myr. We find that the z-velocity 
shows a minimum at ~ 1.5td, when it reaches a distance of 



n{j=1 cm^ , Z{,=500 pc 
n{,=a.1 cm^ , Z{,=500 pc 
no=1 cm^ , Zo=200 pc 
3.1 . zJ^OO pc 




1000 



Figure 3. The ratio min/cs of the 2:— velocity of the top of the 
bubble to the sound speed of the ambient gas at 10"* K is plotted as 
a function of C the mechanical luminosity, and Nqb, the number 
of SNc responsible for the bubble. Different lines correspond to 
different values of mid-plane gas 

number density (1,0.1) cm^"^ and scale heights (200, 500) pc. 



the scale height. We denote this minimum value of z— velocity 
as Uz,min, and refer to this epoch as the 'stalling epoch' in our 
discussion below. 

Figure [3] shows the Mach number at stalling height, as 
a function of C, the mechanical luminosity (which scales 
oc Nob)- Interestingly, superbubbles with Mach number (at 
stalling height) of order less than unity can be triggered by 
even a single SN. These, in principle, can accelerate later and 
therefore breakout of the disk. However, as we shall see later 
with our simulations, there is a minimum number of SNe 
needed for superbubbles to breakout of the disk, particularly 
for high density disks. We also find from Fig |3] that in order 
to achieve a Mach number at stalling height of order 5, one 
needs £ ^ 7 x 10^* erg s"\ for no = 1 cm"^ and zo = 500 
pc. This is larger than the estimate of Koo & McKee (1992), 
and Mac Low & McCray (1988), because of the inclusion of 
radiative loss from the shell. If we consider Vz.min/cs ^ 5 as 
the breakout condition, then we find that larger densities and 
scale heights put more stringent condition on the bubble to 
breakout. 

Next we plot in Figure 3] the minimum Mach number 
as a function of the surface density of Nob, considering the 
surface area of the bubble at the stalling height. Note that we 
are not concerned with the mean surface density of SFR in 
the disk galaxy here. The energy injection considered here is 
localized, but the relevant surface area Eis far as an emerging 
superbubble is concerned, is the area of the bubble in the 
plane of the disk at the point of breaking out. We find that 
for the surface density of energy deposition the analytic curves 
become independent of the scale height and depend only on 
the gas density and number of SNe. This is because the area 
of a superbubble in the plane parallel to the disk, scales with 
Zo, and is a constant for a given scale height. We find that for 
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le-05 0.0001 




de 



(20) 



(21) 



10000 



Figure 4. Mach number of the top of the bubble at the mini- 
mum velocity point is plotted as a function of Nqb divided by the 
cross-sectional area of the bubble at the stalling height, for ana- 
lytical results and for Kompaneets simulations. Analytical results 
are shown for different values of mid-plane gas number densities 
(1,0.1) cm""^ and scale heights (200,500) pc, whereas simulation 
results for Kompaneets runs arc shown for no = 0.1, 1 cm^"^ and 
scale height zq = 200 pc. 



, scale height of 500 pc, the threshold surface density of SNe 



is Nob ~ 1000 kpc" 



5 NUMERICAL SIMULATIONS 

In addition to analytic estimates and approximate calcula- 
tions, we have performed 2-D axisymmetric hydrodynamic 
simulations of breakout using the ZEUS-MP code (Hayes et 
al. 2006). ZEUS-MP is a publicly available, second-order ac- 
curate Eulerian hydrodynamics code. We have carried out 
two sets of simulations; the first set compares numerical sim- 
ulations with the analytic Kompaneets calculation of strong 
shocks in stratified atmospheres (hereafter these runs are re- 
ferred to as 'Kompaneets runs'); and the second set of calcu- 
lations use a more realistic setup, such as disk gravity, mass- 
loading of the ejecta, for shock (superbubble) breakbout in 
starforming galaxies (hereafter these runs will be called 'real- 
istic runs'). 

In this section we introduce the equations that we solve 
numerically, the initial and boundary conditions, and the 
choice of setup parameters. The simulations are run using 
the 2-D axisymmetric, spherical polar (r, 0, (j)) coordinates. 



5.1 Governing Equations 

We solve the following standard Euler's hydrodynamic equa- 
tions including cooling, external gravity, and mass and energy 
loading at inner radii. 

^ = -pV.v+5p(r), (19) 



- = -q-(n,T) + Se(r), 

where d/dt = d/dt + v.V is the Lagrangian derivative, p 
is the mass density, v is the fluid velocity, p is the thermal 
pressure, e = p/(7 — 1) is the internal energy density (we use 
7 = 5/3 valid for an ideal non-relativistic gas), g = —sgn{z)gi 
(sgn[2]=±l for z ^ 0) is the constant external gravity point- 
ing towards the 2 = plane, q~ = neniA{T) is the cool- 
ing term due to radiation where Ue and Ui are the electron 
and ion number densities, A(T) is the cooling function (solid 
line in Fig. 1 and Eq. 12 in Sharma et al. 2010). There are 
source terms in the mass and internal energy equations {Sp, 
Se)- These terms are non-zero and constant only within rin, 
a small injection radius within which supernovae pump mass 
and energy into the interstellar medium. Note that the mass 
loading {Sp) and external gravity (g) terms are used only for 
the realistic simulations and are set to zero for the Kompa- 
neets runs. 

The energy source function is chosen to mimic the en- 
ergy input by supernovae, Se = £/[(47r/3)r-fn], where £ — 
EsnNob /t* = 6.3 X lO^^A^os erg s~^ is the supernova heat- 
ing rate, Esn = 10^^ erg, t» = 50 Myr is the average lifetime 
of main sequence OB stars, and Nob is the number of OB 
stars. The mass-loading source function Sp is chosen as Sp — 
M/[(47r/3)rfn] where M = l3Rf{L/4 x lO^'^erg s-i)M0yr-^ 
where Rf is the return-fraction {— 0.3) and /3 = 3, that in- 
cludes the effect of stellar winds, as inferred by Strickland 
& Heckman (2009) in the case of M82. Tables [T] and [g] list 
the parameters for our Kompaneets and realistic simulations 
respectively. 

We implement energy injection by assuming the de- 
posited energy to be thermalized within a radius rin, which we 
determine from the condition that the corresponding analytic 
solution of superbubble radius enters the Sedov- Taylor phase. 
We assume a mass loading 5 Nob Mq .The superbubble en- 
ters the Sedov- Taylor phase when the ejecta mass equals the 
mass swept up by the shell. We choose this radius to be our 
rin because before this phase, most of the energy of the su- 
perbubble is in kinetic form, and the assumption of most of 
the energy being thermalized is appropriate only in the Sedov- 
Taylor phase. Moreover, rin should be smaller than the radius 
at which the shock becomes radiative. Tables [1] and [2] show 
the values of rin for different simulations. 



5.2 Initial and Boundary conditions 

We have used ZEUS-MP in spherical polar (r, 9, 0) coordi- 
nates. We flx the inner radial boundary (the mass and en- 
ergy injection radius) at rmin "C J"in, and the outer bound- 
ary at Tmax = 3-10 kpc, depending on the distance reached 
by the superbubble in 1.23 x 10^^ s (39.3 Myr) , the maxi- 
mum time for which we run the simulations. For 8 — <j) coor- 
dinates, 6 goes from to tt, and ((> goes from to 27r. We use 
a logarithmically spaced grid in the radial direction such that 
there are equal number of grid points in [rmin, ('"min'"max)'^^^] 
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Table 1. Parameters for Kompaneets runs {C = 6.3 X 10^^ erg 



no (cm ^) 


Nob 


Tin (pc) 


fmin (pc) 


r-max (pc) 


0.1 


1 


10 


5 


2500 


0.1 


10 


21 


10 


2500 


0.1 


100 


44 


30 


2500 


0.1 


300 


63 


40 


3000 


0.1 


1000 


94 


70 


3500 


1.0 


1 


5 


3 


2500 


1.0 


10 


10 


5 


2500 


1.0 


100 


20 


10 


2500 


1.0 


300 


29 


15 


3000 


1.0 


1000 


44 


30 


3000 



and [(rminrmax)^''^, rmax]; the grid is uniformly spaced in the 
other directions. The resolution adopted for our simulation is 
512 X 256 X 1, in the r, 6, (p directions (although we have used a 
higher resolution of 1024 x 1024 x 1 for our study of thermal 
instability in the relevant cases). Outflow boundary condi- 
tions are applied at the outer radial boundary. Inflow-outflow 
boundary condition is applied at the inner radial boundary 
such that mass is allowed to leave or enter the box. Reflective 
boundary conditions are imposed ai 8 = 0, tt, and periodic 
boundary conditions are applied in the (f) direction. 

The initial conditions in Kompaneets and realistic runs 
are different. Both set of runs have an initial temperature 
of 10* K corresponding to the stable WIM. In Kompaneets 
runs the density is stratified in the vertical {z—) direction as 
p{t = 0) oc e~'^^'^", where zq is the scale height. Thus, the 
initial state is not in dynamical equilibrium. However, since 
the sound speed is very small, the evolution occurs because 
of fast energy injection in the center. We have verified that 
the results are the similar as for simulations with a constant 
initial pressure. For the Kompaneets runs, we have used a 
scale height of zo — 200 pc. 

For realistic runs the initial ISM is symmetric with re- 
spect to the vertical direction, with p{t = 0) oc e^'^''^^", where 
the scale height is determined self-consistently for an isother- 
mal gas in hydrostatic equilibrium; i.e., the strength of the 
constant gravitational acceleration is chosen to be g = c^/zo 
where Cs = ^/kTJJjm^ {k is Boltzmann constant, fj, is the 
mean particle mass, and rup is proton mass) is the isothermal 
sound speed and zq is the scale height. 



5.3 Kompaneets runs 

We first describe the results of our Kompaneets runs, of super- 
bubbles in a stratified atmosphere without external gravity or 
mass-loading. Figure [4] shows the variation of the minimum 
Mach number of the top of the superbubble as a function of 
the surface density of energy injection in the disk, for a scale 
height of 200 pc and two values of ambient density, no = 0.1 
and 1 cm~^. We find that the analytical results overestimate 



the Mach number of the superbubbles compared to the simu- 
lations by a factor of order ~ 2, because they do not take into 
account the loss of energy through pdV work. The inclusion 
of radiative loss increases the difference, because the analyt- 
ical estimate of energy loss described in the previous section 
is based on simplified assumptions. 

5.4 Realistic runs 

Next we describe simulations that includes vertical disk grav- 
ity and mass loading. We study the case of ambient gas at 
T = 10* K, with mid-plane densities no =0.1 and 1 cm~'^, 
and scale heights zo = 100 and 500 pc. Our choice of parame- 
ters essentially brackets the possible range of gas density and 
scale height in disk galaxies. For example, the distribution 
of the warm ionized gas in Milky Way has been observed 
to have an exponential profile with no ~ 0.01-0.03 cm~'^ 
and Zo ~ 400-1000 pc (Reynolds 1991; Nordgren et al. 1992; 
Gaensler et al. 2008). Since the dynamics of the superbubbles 
depends on the mid-plane density and scale height, and not 
on the temperature (which only dictates the sound speed and 
therefore the Mach number), we should also consider the ver- 
tical distribution of HI. Dickey & Lockman (1990) found that 
the vertical distribution of HI in Milky Way is best described 
by two Gaussians of central densities 0.4 and 0.1 cm"'', with 
full width half maxima of 210 and 530 pc, and an exponen- 
tial with central density 0.06 cm~^ and a scale height of 400 
pc. The combined distribution has a FWHM of 230 pc and 
a central density of 0.57 cm"''. These values are close to the 
ranges of density and scale heights used in our simulations. 
However, our main goal in this paper is to study the eflect 
of superbubbles in stratified disks in general, and therefore, 
instead of attempting to simulate the Milky Way or any par- 
ticular galaxy, we chose to run a set of idealized simulations 
with exponential envelopes. Furthermore, in a combined dis- 
tribution of a Gaussian and an extended exponential enve- 
lope, the eflect of the former is limited to the central regions, 
whereas the latter aflects the dynamics of the superbubble at 
the point of emerging from the disk and beyond, which is the 
focus of our paper. 

We also use smaller scale heights in our simulations. The 
scale height near the centres of galaxies is smaller than that in 
the outer regions, because of deeper gravitational potentials 
in the central regions. Also Dalcanton, Yoachim, Bernstein 
(2004) found that the HI scale height of disk galaxies varied 
with the rotation speed (or, equivalently, the galactic mass). 
Dwarf spirals with rotation speed 50 km s~^ have zo ~ 200 
pc whereas larger galaxies (with rotation speed in excess of 
120 km s"^) have zo = 500-1000 pc. Also, as Basu et al. 
(1999) have found, the scale height encountered by Milky Way 
superbubbles such as W4 is rather small (^ 100 pc). 

We first find that unlike in the analytical case, where su- 
perbubbles ultimately break out of the disk sooner or later, 
irrespective of the energetics, the realistic simulation runs 
show that for high density disk material (no ^ 1 cm"^) , 
superbubbles keep decelerating for ever for a surface density 
of OB stars ~ 100(zo/100pc) kpc~^, or equivalently a surface 
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Table 2. Parameters for Realistic runs 





77 n ( rm ^ ^ 


^' Urf 


T- (nc) 
' 111 Vl-^^/ 


r • (vic) 
' mm VH^/ 


T"™^,. (nc 
' max Vl-^ 


100 


0.1 


1 


10 


5 


1000 


100 


0.1 


10 


21 


10 


2500 


100 


0.1 


100 


44 


10 


2500 


100 


0.1 


300 


63 


10 


2500 


100 


0.1 


1000 


94 


50 


2500 


100 


1 


100 


20 


10 


2500 


100 


1 


300 


29 


15 


2500 


100 


1 


1000 


44 


30 


2500 


100 


1 


2000 


55 


40 


2500 


100 


1 


3000 


63 


40 


2500 


500 


0.1 


10 


21 


10 


2500 


500 


0.1 


100 


44 


10 


2500 


500 


0.1 


300 


63 


30 


3500 


500 


0.1 


1000 


94 


50 


3500 


500 


0.1 


3000 


135 


110 


3500 


500 


0.1 


10000 


201 


160 


3500 


500 


0.1 


50000 


344 


300 


12000 


500 


0.1 


100000 


433 


400 


12000 


500 


1 


1000 


44 


30 


2500 


500 


1 


2000 


55 


40 


2500 


500 


1 


3000 


63 


40 


2500 


500 


1 


5000 


75 


50 


3500 


500 


1 


10000 


94 


70 


3500 


500 


1 


100000 


201 


150 


5500 



E 




time (Myr) 



Figure 5. Velocity of the topmost point of the bubble is plot- 
ted against time for Nqb = 1000, but for different combina- 
tions of scale height (zq = 100, 500 pc) and mid-plane gas den- 
sity (no = 0.1, 1 cm~^). The horizontal lines in each case shows 
(1/5)(£/po2q)^''^, the expected scaling. 

density of mechanical luminosity of ^ 6 x 10'*'^ (zo/lOO pc) erg 
cm^^ s^^. In other words superbubbles never break out of the 
disk in these cases. The corresponding energy injection sur- 
face density is ~ 2-5 x 10~^ erg cm~^ s~^ For lower density 
ambient gas, no ~ 0.1 cm"'', however, even a single SN event 
can drive a bubble through the disk. We note that this limit is 
consistent with that found by Silich & Tenorio-Tagle (2001) 
for a Milky Way type disk. 

In the case of a superbubble breaking out of the disk, 



there are differences in the way they evolve depending on the 
energy injection rate. We show the evolution of the speed of 
the topmost point of the bubble as a function of time for four 
cases in Figure (5] for two mid-plane densities (no = 0.1, 1 
cm~^), and two scale heights {zq = 100,500 pc), all for a 
surface density of OB stars of 1000 kpc"'^. The curves show 
that the bubbles show acceleration after breakout of the disk 
only for the case of low density and small scale height ( see 
the curve at the top-left corner, for no = O.l cm~^, zq = 100 
pc). In other cases, for disk column density ^ 3 x 10^^ cm~^, 
the bubbles either coast along with the the speed that they 
reach at the breakout, or decelerate to some extent, for a 
considerable period of time before they start accelerating after 
reaching a distance of several scale heights. The curves show 
that the speed at the stalling height, or the minimum speed 
of the bubbles, is an important characteristics of the bubble 
dynamics. It is important because this is the characteristic 
speed with which the bubble sweeps most of the extra-planar 
region of the halo. Also, since the bubble begins to accelerate 
only after reaching a distance of a few times the scale height, 
the corresponding Rayleigh- Taylor instability should not set 
in at the scale height, but at a much larger distance. We shall 
re-visit this point in the next section on instabilities. In some 
cases, the curves show a deceleration at late times. This is due 
to the formation of clumps in the shell from radiative cooling, 
which often sink through the hot gas owing to gravity. 

We have found that typically the minimum speed 
Vz,min ~ (l/5)(£/poZo)^''^ ~ zo/{5td), where td is the dy- 
namical time defined earlier. These values are shown as hori- 
zontal lines in Figure [S] for respective cases. It is easy to see 
that in case of little radiation loss, the speed of the bubble at 
the time of reaching the scale height is ~ {3/5){£,/ pqZq)^^^ , 
as expected from the self-similar evolution of a bubble (r ~ 
{£t^ / poy^^)- Our simulations show that the actual speed is 
roughly a third of this value, and therefore shows the impor- 
tance of radiative loss in the dynamics of superbubbles. As 
analytically derived earlier, radiation losses remove as much 
as half of the total energy of the superbubbles. We recall that 
for an ambient medium with a given temperature, the dimen- 
sionless quantity defined by Mac Low & Norman (1988) is 
D ~ {5vz,min/cs)^ , SO that their condition of D ^ 100 for 
break out corresponds to a minimum Mach number of order 
unity. 

We show the resulting value of minimum Mach number of 
superbubbles for different no and zo in Figure [G] as a function 
of surface density of energy injection. The curves show that in 
terms of energy injection or SNe surface density, the crucial 
parameter is the mid-plane gas density, which separates the 
curves, as was also indicated by our analytical results. Super- 
bubbles with a given surface density of energy injection find 
it easier to break out of disks with lower mid-plane density. 
However, scale height also makes a small difference unlike in 
the analytical calculations; a higher energy density is required 
to clear a thicker disk. 

The important features of our results as shown in Figure 
|6]are: 

• As mentioned above, the condition for a break out from 
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Figure 6. The minimum Mach number of the top of the bubble in 
our reaUstic runs are shown as a function of Nqb per kpc~^, and 
C/irr'^ (erg cm~^ s"'^), for rio = 0.1, 1 cm"'' and zq = 100, 500 pc. 
Note that, for no = 1 cm~^, the shocks stall for a surface density 
of OB stars ^ 500 kpc-^. 



a dense ambient medium with gas density of no = 1 cm~ 
is an energy injection rate surface density of 2-5 x 10^^ erg 
cm~^ s~^. For lower gas densities, the required rate density 
is ~ 10~^ erg cm~^ s~^. The corresponding Mach number for 
these superbubbles that can be as low as of order unity. 

• Superbubbles that can break out with a larger Mach 
number of ~ 3-5 corresponds to ~ 1000 Nqb kpc^^, or an 
energy injection surface density of lO^** erg s^^ cm~'^, for the 
most realistic spiral disks, with no = 0.1 cm~"^, zo = 500 pc, 
or no = 1 cm""^, zo = 100 pc (which for Milky Way case 
describes either the warm extra-planar or cold gas). We note 
that the largest OB associations have ~ 10* Mq (McKee 
& WiUiams 1997), and with a cross-sectional area of order 
7r(7r2o)^ (since rmax ~ t^zq asymptotically; see equation [9]) , a 
superbubble blown by such a large OB association can only 
have ^ 10^ SNe kpc"^ for a Salpeter IMF. Therefore we 
can conclude that (a) only the largest of the OB associations 
can produce a bubble that can break out of Milky Way-type 
disks, (b) this also corresponds to the minimum energy injec- 
tion rate of 10~* erg s^'^ cm~^ as observed by Dahlem et al. 
(1995) for the existence of radio emitting halo gas, and (c) 
larger ISM density or scale height would require more than 
one OB association to produce a superbubble or adjacent mul- 
tiple bubbles that can coalesce and grow together. Recent 
simulations show that cosmic rays can stream through ISM 
gas to considerable heights above the disk, and break out of 
superbubbles can provide such channels (Uhlig et al. (2012). 

• As explained earlier, the minimum speed of 3-5 Cs ~ 30 
km s~^, for an ambient gas at 10* K, also corresponds to the 
case where the hot (and multiphase gas; see next section on 
instabilities and gas cooling) interior gas can sweep up to a 
height of ~ 1 kpc within a time period of ~ 50 Myr, the time 
scale over which OB stars explode and keep injecting energy 



in the bubble. Combined with the result mentioned above, 
we can conclude that an energy injection rate of 10"* erg s~* 
cm^^, or ~ 1000 Nqb kpc^^ not only can produce a bubble 
that can break out of the disk but also fill the halo up to a 
height of order ~ 1 kpc. 

• If we insist on a larger Mach number at stalling height, 
to be 5-10, then the energy injection rate becomes ~ 10"'^ 
erg s^* cm^^, with ~ 2x 10* Nqb kpc^^. Using a time scale of 
~ 50 Myr of OB stars, the corresponding SFR surface density 
for a Salpeter IMF is ~ 0.06 M© yr~* kpc"'^. If superbubbles 
seed galactic outflows, then the gas speed is required to be a 
few hundred km s~^, and the Mach number at stalling height 
is needed to be much larger than ten, and the corresponding 
requirement on SFR surface density increasing to ~ 0.1 Mq 
yr^'^ kpc~^, the observed threshold. Therefore, the Heckman 
(2000) threshold ( ~ 0.1 Mq yr^^ kpc^^) for superwinds cor- 
responds to a larger requirement on the part of superbubbles, 
of not only breaking out of disks but doing so with a large 
Mach number. 



6 THERMAL AND RAYLEIGH- TAYLOR 
INSTABILITY 

The focus till now was on the important Vz/cs parameter 
(the minimum Mach number of the shell) which determines 
the fate of the superbubble after it crosses the scale height. In 
this section we discuss the role of different instabilities, in par- 
ticular Ray leigh- Taylor and thermal instabilities, in our 2-D 
breakout simulations. When the superbubble reaches about a 
scale height, the shock is generally believed to accelerate ow- 
ing to the decrease in pressure. This should lead to the onset 
of the Rayleigh- Taylor (RT) instability, as has been invoked in 
previous analytical works (e.g., Koo & McKee 1992) and seen 
in numerical simulations (e.g., Mac Low, McCray & Norman 
1989). However, as mentioned earlier, our simulations show 
that superbubbles do not accelerate until after they reach a 
distance of several scale heights ( as was also found by Fer- 
rara & Tolstoy 2000). Therefore RT instability occurs at a 
distance much larger than the scale height. Also we find that 
before the onset of RT instability, the superbubble expand- 
ing in the disk suffers from thermal instability in the early 
stages of its evolution. This instability leads to clumping and 
fragmentation of the shell of the superbubble well in advance 
of the RT instability, and can therefore affect the outcome of 
the RT instability. 

Figure [7] shows the 2-D snapshots of temperature at two 
different times for our fiducial high resolution run {Nqb = 
5000, no = 1 cm~^, zo = 500 pc). Figure [6] indicates that the 
minimum Mach number for this case is ~ 2 and the bubble is 
just about able to break out within the starburst timescale. 
The temperature snapshot at early time (9 Myr), when the 
bubble has just reached the scale-height, shows that the bub- 
ble is roughly spherical. The radiative shell seems to develop 
corrugations where the hot bubble gas and the radiatively 
cooled shocked gas interpenetrate. The shell is at ~ 10* K 
(the same as the ambient ISM temperature), the tempera- 
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Figure 7. Temperature contours (colour coded) for a superbubblc 
with Nob = 5000, no = 1 cm~'^, zq = 500 pc, at t = 9 Myr, when 
the top of the bubble has reached a distance of the scale height 
(left panel), at 39.3 Myr, when it has reached a distance ~ 3zo 
(middle panel). The rightmost panel shows the case of the same 
superbubble without radiative cooling at t = 39.3 Myr, the same 
evolutionary epoch as the middle panel. 



ture below which the cooling function drops suddenly and 
the gas becomes thermally stable. The dense shell is more 
clearly seen in the density snapshots of Figure [S] The corru- 
gations are definitely driven by radiative cooling because the 
run without radiative cooling shows a smooth shell (the third 
panel in Figs. [7] and 

While the fragments of cold shell are confined to the bub- 
ble boundary at early times, the cold gas lags behind the hot 
gas at later times because the hot gas is pushed out by super- 
nova heating. The cold blobs are only pushed out because of 
the drag force due to the hot gas but eventually trail behind. 
The cold blobs embedded in the hot gas are reminiscent of the 
cold multiphase filaments observed in galactic outfiows, such 
as M82. Since in our simulations cold gas leaves the simulation 
box from the inner boundary, all the cold blobs embedded in 
the hot bubble come from the fragmenting cold shell. In real- 
ity, some cold gas from the cold star-forming regions can also 
be uplifted by the hot gas. At late times, in the runs with 
cooling, there are some signs of bubble breaking out because 
of RT instability close to the polar regions. All such signa- 
tures of RT instability are missing in the run without cooling 
(panel 3). This is mainly because RT instability in the run 
with cooling is seeded with large amplitude perturbations by 
corrugations caused by shell cooling. 

In order to assess the relative importance of thermal and 
RT instabilities, in Figure [9] we show the free-fall time (iff = 
\j2zl(vz + <?), where z and are the height and acceleration 
of the shell, and g is the acceleration due to gravity; and the 
cooling time (tcooi = l.ZkT /nA) of the shell as a function of 
time for runs corresponding to Figures [7] & [S] We use the 
position of the outermost densest part to identify the shell 
position. In the left panel of Figure [9l we show the case of 
Nob = 5000 no — 1 cm"'^, zq = 500 pc. We expect the shell 
to cool radiatively if fcooi is shorter than time. And indeed, 
the radiative cooling time is shorter than time at early times. 
This is consistent with the cooling and fragmentation of the 
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Figure 8. Density contours for the same cases as in Fig [7] Here, 
fragmentation of the shell is clearly seen in the run with cooling. 
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Figure 9. The free-fall and cooling timescalcs for the shell mate- 
rial are plotted against time, for two examples with Nqb = 5000, 
20 = 500 pc, and no = 1 cm~'^ (left panel), no = 0.1 cm~^ (mid- 
dle panel). The grey lines show the time elapsed in each cases for 
comparison. The right panel shows the case of no radiation cooling 
for no = 1 cm^^. The leftmost and rightmost panels correspond 
to the runs shown in Figs Ffl and l8l 



dense shell seen in Figs.[7]&[Sl One point of caution: we should 
ideally plot the cooling time of the shell assuming the shell 
temperature and density corresponding to an adiabatic shock 
because cooling will happen if this timescale is short. Here 
we are plotting the cooling time of the shell, which for the 
left panel case, has already cooled to low temperatures. Since 
cooling time increases sharply below 10* K, tcooi is barely 
smaller than time in the left panel of Figure |9l At later times 
fcooi becomes longer than time and we do not expect the 
newly accumulated shell material to cool. The RT timescale 
(« ts) is always longer than time for the fiducial run. The 
free-fall time increases initially as the shock slows down until 
a scale height. After that the shock moves at a small Mach 
number ~ 2. This is consistent with the fact that we do not 
see vigorous RT instability in Figures [7| & [S] 

The middle panel of the Figure[5]shows various timescales 
for a midplane density of no =0.1 cm^'^ The cooling time for 
this case is shorter than the cooling time for the higher den- 
sity case. This seems inconceivable given the higher density 
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and efficient cooling for the run in the left panel. This dis- 
crepancy arises because although the density for the no = 0.1 
cm~^ is smaller, the temperature of the post-shock gas is 10^ 
K, where the cooling function peaks. Consequently the cool- 
ing time is shorter than the higher density run. For compari- 
son, we have also plotted the cooling and free-fall timescales 
for the runs without cooling in the right panel. The density 
and temperature snapshots for this run do not show cooling- 
induced fragmentation. 

We have shown in Figure [6] by darkened points the cases 
in which tcooi is always less than tfi, for different values of 
C/nr'^, zo and no- We find that these cases mostly appear 
for which, roughly, 10 ^ ^'^.min 3, except for the case of 
Zo = 500 pc and no = 1 cm^^, for which there is a cross- 
over point in time after tcooi ^ iff. We note that this range of 
Vz,min corresponds to a case in which the shell temperature 
(Ts) remains in the range of 2 x 10"* ^ Ts ^ 10®, where 
the cooling function peaks. This implies a range in Nob for 
which thermal instability is imporant. In the low Nqb limit, 
the shock is not strong enough and Ts ^ 10'' K, and in the 
high Nob case, the shock is very strong (Ts > 10® K) and tg 
(RT timescale) is shorter than tcooi at late times. 

We are therefore led to conclude that superbubbles are 
affected not only by RT instability but also by thermal insta- 
bility, depending on the density and energy injection. This im- 
plies that the fragmentation of the bubble shell that releases 
the hot interior gas into the halo occurs under the combined 
effects of thermal instability at early times and RT instabil- 
ity at late times if the Mach number at stalling epoch is large 
enough. 



7 DISCUSSION & SUMMARY 

Superbubbles with fragmented shells are believed to ulti- 
mately form 'chimneys' (Norman & Ikeuchi (1989), which 
connect the halo gas to the processes in the disk in different 
ways. Apart from transporting hot gas to the halo, chimneys 
provide a natural channel for Lyman continuum photons from 
hot stars in the disk to reach the diffuse ionized medium of 
the Reynolds layer (Reynolds 1991; Dove & ShuU 1994). It 
is however important for the superbubble shells to fragment 
before the main sequence life times of O stars for a substan- 
tial fraction of ionizing radiation to escape the disk (Dove, 
ShuU, Ferrara 2000). This implies a fragmentation time scale 
of 3-5 Myr, which is comparable to the dynamical timescale 
{td ~ zf/'^^ipo/Cy^'^), for superbubbles with C ~ 10^** erg 
(corresponding to Nob ~ 200), typical disk parameters. This 
is the energy scale for the largest of the OB associations, and 
as our results show superbubbles with smaller energetics find 
it hard to pierce through the disk, unless the OB association 
is located much above the mid-plane level. 

In other words, for superbubbles to act as effective con- 
duits of ionizing radiation for the halo, or for the intergalactic 
medium (at high redshift, in the context of the epoch of reion- 
ization), the superbubbles need to fragment roughly around 
the time when they reach a scale height. This is unlikely to 



happen only through RT instability as superbubbles do not 
accelerate until reaching a distance of several scale heights. 
Also, as de Avillez & Breitschwerdt (2005) have discussed 
on the basis of simulations of a magnetized ISM, superbub- 
ble shells can stabilize against RT instability in the presence 
of magnetic fields. In this regard, the clumping of the shell 
from thermal instability at an early phase of evolution of the 
superbubble can be important. 

We have studied the evolution of superbubbles in strati- 
fied disks analytically and with simulations. Our results can 
be summarised as follows: 

• Our analytic calculations show that radiation losses are 
important for superbubble dynamics. Radiation loss is more 
important for superbubbles with continuous energy injection 
than a supernova remnant of similar total energy. We estimate 
a fraction 4/9 of the total energy being being radiated away. 
We have further checked our analytical results, based on a 
modified Kompaneets approximation that takes into account 
radiative cooling, with numerical simluations. We found that 
analytic results differ from that of idealized simulations by 
a factor ~ 2. The results obtained by the analytical means 
therefore provide a useful benchmark to compare with realis- 
tic simualtions. Also, for disks with large gas density, with 
no ^ 1 cm^'^, superbubble breakouts are not possible for 
surface density of OB stars ^ 100(zo/100 pc) kpc~'^, or an 
equivalent energy injection surface density of ^ (2-5) x 10~^ 
erg cm~^ s~'. 

• Superbubbles that emerge from the disk with Mach num- 
ber of order 2-3 require an energy injection rate of ~ lO"'' 
erg cm~^ s~', corresponding to explosions triggered by the 
largest OB associations with 10* Mq. This energy injection 
scale corresponds to disk galaxies with synchrotron emitting 
gas in the extra-planar regions. 

• Vigorous superbubbles that break out of the disk with 
sufficiently large Mach number (^ 10) , correspond to an en- 
ergy injection rate of ~ 10"'^ erg cm^'^ s^', or equivalently, 
a SFR surface density of ~ 0.1 M© yr^' kpc^'^. These super- 
bubbles require more than one OB associations to produce 
and sustain their dynamics, and this energy injection scale 
corresponds to (a) the existence of multiphase gas in the halo 
of disk galaxies, and (b) the Heckman threshold for the onset 
of superwinds. 

• Superbubbles do not accelerate until reaching a vertical 
distance of a few scale heights (of order ~ 2), which implies 
that RT instability helps to fragment the shells not at a dis- 
tance of a scale height but at a much larger height. Also, we 
find that for typical disk parameters, thermal instability acts 
on the shell at the early stages of superbubble evolution, and 
forms clumps and fragments in the shell, much before the shell 
is acted upon by RT instability. Radiative cooling therefore 
manifests in seeding thermal instability, which has important 
implications for the clumping of superbubble shell and pro- 
ducing channels of leakage for ultraviolet radiation into the 
halo. 
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